(from http://journals.aps.org/prstab/abstract/10.1103/PhysRevSTAB.17.014001)
In [1]:
!date
%matplotlib inline
In [21]:
from __future__ import division
def fareySequence(N, k=1):
"""
Generate Farey sequence of order N, less than 1/k
"""
# assert type(N) == int, "Order (N) must be an integer"
a, b = 0, 1
c, d = 1, N
seq = [(a,b)]
while c/d <= 1/k:
seq.append((c,d))
tmp = int(math.floor((N+b)/d))
a, b, c, d = c, d, tmp*c-a, tmp*d-b
return seq
def resonanceSequence(N, k):
"""
Compute resonance sequence
Arguments:
- N (int): Order
- k (int): denominator of the farey frequency resonances are attached to
"""
a, b = 0, 1
c, d = k, N-k
seq = [(a,b)]
while d >= 0:
seq.append((c,d))
tmp = int(math.floor((N+b+a)/(d+c)))
a, b, c, d = c, d, tmp*c-a, tmp*d-b
return seq
def plotResonanceDiagram(N, figsize=(10,10)):
import matplotlib.pyplot as plt
ALPHA = 0.5/N
plt.figure(figsize=figsize)
ticks = set([])
for h, k in fareySequence(N, 1):
ticks.add((h,k))
for a, b in resonanceSequence(N, k):
if b == 0:
x = np.array([h/k, h/k])
y = np.array([0, 1])
elif a== 0:
x = np.array([0, 1])
y = np.array([h/k, h/k])
else:
m = a/b
cp, cm = m*h/k, -m*h/k
x = np.array([0, h/k, 1])
y = np.array([cp, 0, cm+m])
plt.plot( x, y, 'b', alpha=ALPHA) # seqs. attached to horizontal axis
plt.plot( y, x, 'b', alpha=ALPHA) # seqs. attached to vertical axis
# also draw symetrical lines, to be fair (otherwise lines in the
# lower left traingle will be duplicated, but no the others)
plt.plot( x, 1-y, 'b', alpha=ALPHA)
plt.plot(1-y, x, 'b', alpha=ALPHA)
plt.xlim(0, 1)
plt.ylim(0, 1)
plt.xticks([h/k for h,k in ticks], [r"$\frac{{{:d}}}{{{:d}}}$".format(h,k) for h,k in ticks], fontsize=15)
plt.yticks([h/k for h,k in ticks], [r"${:d}/{:d}$".format(h,k) for h,k in ticks])
plt.title("N = {:d}".format(N))
Try it!
Note: in the original paper there was a minor mistake. Eq. (8) read
$$ \Bigg( \Big\lfloor \frac{N+b+a}{d} \Big\rfloor c - a , \Big\lfloor \frac{N+b+a}{d} \Big\rfloor d - b \Bigg) $$but it should read
$$ \Bigg( \Big\lfloor \frac{N+b+a}{d+c} \Big\rfloor c - a , \Big\lfloor \frac{N+b+a}{d+c} \Big\rfloor d - b \Bigg) $$I've contacted Rogelio Tomás (the author) and he agreed with the correction (an erratum will be sent to the publication).
In [22]:
N = 5
for k in set([k for _,k in fareySequence(N, 1)]):
print "N={}, k={}:".format(N, k)
print "\t", resonanceSequence(N, k)
In [23]:
# from matplotlib2tikz import save as save_tikz
plotResonanceDiagram(5, figsize=(12,12))
# save_tikz('resonanceDiagram_N7.tikz')
def plotSolution(a,b,c,color='r'):
x = [c/a, 0, (c-b)/a, 1]
y = [0, c/b, 1, (c-a)/b]
plt.plot(x, y, color=color, alpha=0.5, linewidth=4)
# plot some example solutions
if True:
# solutions for (x,y) = (0.5, 1)
plotSolution( 4, -1, 1)
plotSolution(-2, 2, 1)
plotSolution(-2, 3, 2)
plotSolution( 4, 1, 3)
plotSolution( 2, 1, 2)
# solutions for (x,y) = (0.5, 0.5)
plotSolution( 3, -1, 1, 'g')
plotSolution(-1, 3, 1, 'g')
plotSolution( 3, 1, 2, 'g')
plotSolution( 1, 3, 2, 'g')
plotSolution( 1, 1, 1, 'g')
In [ ]: